#ifndef _NOR_ORDER_
#define _NOR_ORDER_

#include "../DF_Common/basis_state.h"
#include "../DF_Common/coupling_coeff.h"
#include "../DF_Local/package.h"
#include "../DF_Common/read_data.h" 
#include "../DF_Common/time_count.h"
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
using std::string;

class Norm_Ord {
public:
  Norm_Ord(SingleP &SP_x, TripleP &TP_x, SMatrix_Ch &vec_Tz_x,
           SMatrix_Ch &vec_flag_Tz_x)
      : SP(SP_x), TP(TP_x), vec_Tz_S(vec_Tz_x), vec_flag_Tz_S(vec_flag_Tz_x) {}

  SingleP &SP;
  TripleP &TP;
  SMatrix_Ch &vec_Tz_S, &vec_flag_Tz_S;
  // Package &pack;
  vector<vector<vector<double>>> vec_Tz;
  vector<vector<vector<int>>> vec_flag_Tz;

  int Num_valence_p, Num_valence_n;

  vector<vector<State_DP>> state_T2P_Tz; // [ch][i]
  vector<Channel_DP> channel_T2P_Tz;
  vector<vector<vector<double>>> vec_2f3, vec_2f3_C,
      vec_2f3_P; //[ch_J][i][j] two body matrix element from V3N
  vector<double> vec_1f3, vec_1f3_C, vec_1f3_P,
      vec_1f3_1C1P; // [i] one body matrix element from V3N
  double val_0f3;   // zero body energy from V3N

  vector<vector<double>> rho_vec; //[Tz:0,1,2][p: SP_Tz]

  vector<double> rho_p_vec; // tz = -1, 0, 1.  nn, pn, pp
  vector<double> rho_n_vec; // tz = -1, 0, 1.  nn, pn, pp

  vector<double> rho_p_vec_sp; // tz = -1, 1. n p
  vector<double> rho_n_vec_sp; // tz = -1, 1. n p

  double rho_p_H0; // tz = -1, 1. n p
  double rho_n_H0; // tz = -1, 1. n p

  void build_T2P_Tz(int a_min, int a_max, int b_min, int b_max);
  int FindC_T2P_Tz(int J2ab, int T2ab_z);
  int Find_T2P_Tz(int J2ab, int T2ab_z, int a, int b);

  void cal_val0f3_Tz(int h_core, int h_P_beg, int h_P_end, int Num_p,
                     int Num_n);

  double cal1f3_Tz_C(State_SP &sp_t, int h_core);
  double cal1f3_Tz_P(State_SP &sp_t, int h_P_beg, int h_P_end,
                     vector<double> rho_p_t, vector<double> rho_n_t);
  double cal1f3_Tz_1C1P(State_SP &sp_t, int h_core, int h_P_beg, int h_P_end,
                        vector<double> rho_p_t, vector<double> rho_n_t);

  double cal2f3_Tz_C(State_DP &T2P_i, State_DP &T2P_f, int h_core);
  double cal2f3_Tz_P(State_DP &T2P_i, State_DP &T2P_f, int h_P_beg, int h_P_end,
                     vector<double> rho_p, vector<double> rho_n);
  double cal2f3_Tz_P2(State_DP &T2P_i, State_DP &T2P_f, int h_P_beg,
                      int h_P_end, vector<vector<double>> rho_val);

  void build_vec2f3_Tz_C(int h_core);
  void build_vec2f3_Tz_P(int h_P_beg, int h_P_end, vector<double> rho_p_t,
                         vector<double> rho_n_t);
  void build_vec2f3_Tz_P2(int h_P_beg, int h_P_end,
                          vector<vector<double>> rho_val);
  void build_vec2f3_Tz(int h_core, int h_P_beg, int h_P_end, int Num_p,
                       int Num_n);
  void build_vec2f3_Tz2(int h_core, int h_P_beg, int h_P_end, int Num_p,
                        int Num_n, int flag_P);
  void build_vec2f3_Tz3(int sp_core_Tz);

  void build_vec1f3_Tz_C(int h_core, int p_beg, int p_end);
  void build_vec1f3_Tz_P(int p_beg, int p_end, int h_P_beg, int h_P_end,
                         vector<double> rho_p, vector<double> rho_n);
  void build_vec1f3_Tz_1C1P(int h_core, int p_beg, int p_end, int h_P_beg,
                            int h_P_end, vector<double> rho_p,
                            vector<double> rho_n);
  void build_vec1f3_Tz(int h_core, int p_beg, int p_end, int h_P_beg,
                       int h_P_end, int Num_p, int Num_n);
  void build_vec1f3_Tz2(int h_core, int p_beg, int p_end, int h_P_beg,
                        int h_P_end, int Num_p, int Num_n, int flag_P);
  void build_vec1f3_Tz3(int sp_core_Tz, int sp_P_beg_Tz, int sp_P_end_Tz);

  void print_T2P_Tz(std::string filename);
  void print_channel_Tz(std::string filename);

  void print_vec2f3_Tz(std::string filename);
  void print_vec1f3_Tz(std::string filename);
  void print_norm_ME_C(std::string filename, int h_core, int p_beg, int p_end,
                       double hw, int reg_pow, int Lamb_Cut,
                       vector<int> order_change);
  void print_norm_ME_CP(std::string filename, int h_core, int p_beg, int p_end);
  void print_norm_ME_CP2(std::string filename, int h_core, int p_beg,
                         int p_end);
  void print_with_SecOrd(std::string filename, int h_core_Tz, int p_beg_Tz,
                         int p_end_Tz, vector<double> vec_sec,
                         vector<vector<double>> vec_sec_sub, double hw,
                         int reg_pow, int Lamb_Cut, int E_cut);

  void print_norm_ME_C_Tz(std::string filename, int h_core, int p_beg,
                          int p_end);
  void print_norm_ME_CP_Tz(std::string filename, int h_core, int p_beg,
                           int p_end);

  void checkVsTok_pack_Tz(std::string file1, std::string file2);
  void checkVsTok_norm2f3_Tz(std::string file1, std::string file2);
  void checkVsTok_norm2f3_Tz();

private:
  //
  int Tri(int &j2a, int &j2b, int &j2c);
};

class Norm_Ord_Sec {
public:
  Norm_Ord_Sec(SingleP &SP_x, DoubleP &DP_x, TripleP &TP_x, Read_Data_II &rd_x,
               SMatrix_Ch &vec_Tz_x, SMatrix_Ch &vec_flag_Tz_x)
      : SP(SP_x), DP(DP_x), TP(TP_x), RD(rd_x), vec_3V_Tz_S(vec_Tz_x),
        vec_3V_flag_Tz_S(vec_flag_Tz_x) {}

  SingleP &SP;
  DoubleP &DP;
  TripleP &TP;
  Read_Data_II &RD;
  // Package &pack;
  SMatrix_Ch &vec_3V_Tz_S, &vec_3V_flag_Tz_S;
  // vector<vector<vector<double>>> vec_3V_Tz;
  // vector<vector<vector<int>>> vec_3V_flag_Tz;

  double hbar_omega;

  int E_cut;
  int sp_max;
  int sp_core;
  int sp_P_beg;
  int sp_P_end;
  int sp_Q_beg;
  int sp_Q_end;

  int sp_max_Tz;
  int sp_core_Tz;
  int sp_P_beg_Tz;
  int sp_P_end_Tz;
  int sp_Q_beg_Tz;
  int sp_Q_end_Tz;

  vector<int> sp2rd; //  [SP state index] = Read_Data.J_basis[i].index

  vector<vector<double>> vec_U, vec_UT;
  vector<vector<int>> vec_U_flag, vec_UT_flag;

  vector<vector<vector<double>>> vec_2V2N, vec_2V3N,vec_2V2N_unnorm;

  vector<vector<vector<int>>> vec_2V2N_flag, vec_2V3N_flag;

  vector<double> vec_A1, vec_A2, vec_B1, vec_B2, vec_C1, vec_C2, vec_D1, vec_D2,
      vec_E1, vec_E2, vec_E3, vec_F1, vec_F2;
  vector<double> vec;
  vector<vector<double>> vec_sub;

  void init(double hw, int E_cut, int sp_max, int sp_core, int sp_P_beg,
            int sp_P_end, int sp_Q_beg, int sp_Q_end);

  double cal1f2_Tz_C(int a, int b);
  double cal1f3_Tz_C(int a, int b);
  double cal2f3_Tz_C(State_DP &T2P_i, State_DP &T2P_f);
  void build_U();
  void build_2V();
  void build_2V_Luigi();
  void print_2V_Luigi(string filename);

  double cal_A1(int a);
  double cal_A2(int a);
  double cal_B1(int a);
  double cal_B2(int a);
  double cal_C1(int a);
  double cal_C2(int a);
  double cal_D1(int a);
  double cal_D2(int a);
  double cal_E1(int a);
  double cal_E2(int a);
  double cal_E3(int a);
  double cal_F1(int a);
  double cal_F2(int a);
  void build_sec();

  void build_vec1f3_Tz_C(int h_core, int p_beg, int p_end);

  void build_vec1f3_Tz(int h_core, int p_beg, int p_end, int h_P_beg,
                       int h_P_end, int Num_p, int Num_n);
  void build_vec1f3_Tz2(int h_core, int p_beg, int p_end, int h_P_beg,
                        int h_P_end, int Num_p, int Num_n);
  int pauli_exc_Tz(int a, int b, int c);
  void print_U(string filename1);
  void print_U_all(string filename1);
  void print_2V(string filename1);
  void print_2V_all(string filename1);

private:
  //
  int Tri(int &j2a, int &j2b, int &j2c);
};

class Norm_Ord_Core {
public:
  Norm_Ord_Core(SingleP &SP_x, DoubleP &DP_x, TripleP &TP_x,
                SMatrix_Ch &vec_Tz_x, SMatrix_Ch &vec_flag_Tz_x)
      : SP(SP_x), DP(DP_x), TP(TP_x), vec_3V_Tz_S(vec_Tz_x),
        vec_3V_flag_Tz_S(vec_flag_Tz_x) {}

  SingleP &SP;
  DoubleP &DP;
  TripleP &TP;
  // Package &pack;
  SMatrix_Ch &vec_3V_Tz_S, &vec_3V_flag_Tz_S;
  // vector<vector<vector<double>>> vec_3V_Tz;
  // vector<vector<vector<int>>> vec_3V_flag_Tz;

  double hbar_omega;

  int E_cut;
  int sp_max;
  int sp_core;
  int sp_P_beg;
  int sp_P_end;
  int sp_Q_beg;
  int sp_Q_end;

  int sp_max_Tz;
  int sp_core_Tz;
  int sp_P_beg_Tz;
  int sp_P_end_Tz;
  int sp_Q_beg_Tz;
  int sp_Q_end_Tz;

  vector<vector<vector<double>>> vec_2V3N;
  vector<vector<vector<int>>> vec_2V3N_flag;

  void init(double hw, int sp_max, int sp_core, int sp_P_beg,
            int sp_P_end, int sp_Q_beg, int sp_Q_end);

  double cal2f3_Tz(State_DP &T2P_i, State_DP &T2P_f);

  void build_P();
  void build_Core();

  int pauli_exc_Tz(int a, int b, int c);
  void print_norm_core(std::string filename, int h_core_Tz, int p_beg_Tz,
                       int p_end_Tz, double hw, int reg_pow, int Lambda_Cut,
                       vector<int> oc_o2z);

private:
  //
  int Tri(int &j2a, int &j2b, int &j2c);
};

#endif
